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ABSTRACT 

This  report  considers  the  derivation  of  the  mathematical  model  for  a  missile  autopilot  in  state 
space  form.  The  basic  equations  defining  the  airframe  dynamics  are  non-linear,  however, 
since  the  non-linearities  are  "structured"  (in  the  sense  that  the  states  are  of  quadratic  form)  a 
novel  approach  of  expressing  this  non-linear  dynamics  in  state  space  form  is  given.  This 
should  provide  a  useful  way  to  implement  the  equations  in  a  computer  simulation  program 
and  possibly  for  future  application  of  non-linear  analysis  and  synthesis  techniques, 
particularly  for  autopilot  design  of  missiles  executing  high  g-manoeuvres. 

This  report  also  considers  a  locally  linearised  state  space  model  that  lends  itself  to  better 
known  linear  techniques  of  the  modem  control  theory.  A  coupled  multi-input  multi-output 
(MIMO)  model  is  derived  suitable  for  both  the  application  of  the  modem  control  techniques 
as  well  as  the  classical  time-domain  and  frequency  domain  techniques.  This  is  validated  by 
comparing  the  model  with  the  other  published  results,  and  through  both  open  and  closed- 
loop  systems  simulations.  The  models  developed  are  useful  for  further  research  on  precision 
optimum  guidance  and  control.  It  is  hoped  that  the  model  will  provide  more  accurate 
presentations  of  missile  autopilot  dynamics  and  will  be  used  for  adaptive  and  integrated 
guidance  &  control  of  agile  missiles. 
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Mathematical  Models  for  a  Missile  Autopilot 

Design 


Executive  Summary 

Requirements  for  next  generation  guided  weapons,  particularly  with  respect  to  their 
capability  to  engage  high  speed,  highly  agile  targets  and  achieve  precision  end-game 
trajectory,  has  prompted  a  revision  of  the  way  in  which  the  guidance  and  autopilot 
design  is  undertaken.  This  report  considers  the  derivation  of  the  mathematical  models 
for  a  missile  autopilot  in  state  space  form.  The  basic  equations  defining  the  airframe 
dynamics  are  non-linear,  however,  since  the  non-linearities  are  "structured"  (in  the 
sense  that  the  states  are  of  quadratic  form)  a  novel  approach  of  expressing  this  non¬ 
linear  dynamics  in  state  space  form  is  given.  This  should  provide  a  useful  way  to 
implement  the  equations  in  a  computer  simulation  program  and  possibly  for  future 
application  of  non-linear  analysis  and  synthesis  techniques. 

This  report  also  considers  a  locally  linearised  state  space  model  that  lends  itself  to 
better  known  linear  techniques  of  the  modem  control  theory.  A  coupled  multi-input 
multi-output  (MEMO)  model  is  derived  suitable  for  both  the  application  of  the  modem 
control  techniques  as  well  as  the  classical  time-domain  and  frequency  domain 
techniques.  This  is  validated  by  comparing  the  model  with  the  other  published  results, 
and  through  both  open  and  closed-loop  systems  simulations.  The  models  developed 
are  useful  for  further  research  on  precision  optimum  guidance  and  control.  It  is  hoped 
that  the  model  will  provide  more  accurate  presentations  of  missile  auto-pilot  dynamics 
and  will  be  used  for  adaptive  and  integrated  guidance  &  control  of  agile  missiles. 
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1.  Introduction 


Requirements  for  next  generation  guided  weapons,  particularly  with  respect  to  their 
capability  to  engage  high  speed,  highly  agile  targets  and  achieve  precision  end-game 
trajectory,  has  prompted  a  revision  of  the  way  in  which  the  guidance  and  autopilot 
design  is  undertaken.  Integrating  the  guidance  and  control  function  is  a  synthesis 
approach  that  is  being  pursued  as  it  allows  the  optimisation  of  the  overall  system 
performance.  This  approach  requires  a  more  complete  representation  of  the  airframe 
dynamics  and  the  guidance  system.  The  use  of  state  space  model  allows  the  application 
of  modern  control  techniques  such  as  the  optimal  control  and  parameter  estimation 
techniques  to  be  utilised.  In  this  report  we  derive  the  autopilot  model  that  will  serve  as 
a  basis  for  an  adaptive  autopilot  design  and  allow  further  extension  of  this  to 
integrated  guidance  and  control  system  design. 

Over  the  years  a  number  of  authors  [1-3,  6]  have  considered  modelling,  analysis  and 
design  of  autopilots  for  atmospheric  flight  vehicles  including  guided  missiles.  In  the 
majority  of  the  published  work  on  autopilot  analysis  and  design,  locally  linearised 
versions  of  the  model  with  decoupled  airframe  dynamics  has  been  considered.  This 
latter  simplification  arises  out  of  the  assumption  that  the  airframe  and  its  mass 
distribution  are  symmetrical  about  the  body  axes,  and  that  the  yaw,  pitch  and  roll 
motion  about  the  equilibrium  state  remain  "small".  As  a  result,  most  of  the  autopilot 
analysis  and  design  techniques,  considered  in  open  literature,  use  classical  control 
approach,  such  as:  single  input/  single  output  transfer-functions  characterisation  of  the 
system  dynamics  and  Bode,  Nyquist,  root-locus  and  transient  response  analysis  and 
synthesis  techniques  [5,7].  These  techniques  are  valid  for  a  limited  set  of  flight  regimes 
and  their  extension  to  cover  a  wider  set  of  flight  regimes  and  airframe  configurations 
requires  autopilot  gain  and  compensation  network  switching. 

With  the  advent  of  fast  processors  it  is  now  possible  to  take  a  more  integrated 
approach  to  autopilot  design.  Modern  optimal  control  techniques  allow  the  designer  to 
consider  autopilots  with  high-order  dynamics  (large  number  of  states)  with  multiple 
inputs/ outputs  and  to  synthesise  controllers  such  that  the  error  between  the 
demanded  and  the  achieved  output  is  minimised.  Moreover,  with  real-time 
mechanisation  any  changes  in  the  airframe  aerodynamics  can  be  identified  (parameter 
estimation)  and  compensated  for  by  adaptively  varying  the  optimum  control  gain 
matrix.  This  approach  should  lead  to  missile  systems  that  are  able  to  execute  high  g- 
manoeuvres  (required  by  modem  guided  weapons),  adaptively  adjust  control 
parameters  (to  cater  for  widely  varying  flight  profiles)  as  well  as  account  for  non- 
symmetric  airframe  and  mass  distributions. 

Typically,  for  a  missile  autopilot,  the  input  is  the  demanded  control  surface  deflection 
and  outputs  are  the  achieved  airframe  (lateral)  accelerations  and  body  rates  measured 
about  the  body  axes.  Other  input/ output  variables  (such  as:  the  flight  path  angle  and 
angle  rate  or  the  body  angles)  can  be  derived  directly  from  lateral  accelerations  and 
body  rates. 
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This  report  considers  the  derivation  of  the  mathematical  model  for  a  missile  autopilot 
in  state  space  form.  The  basic  equations  defining  the  airframe  dynamics  are  non-linear, 
however,  since  the  non-linearities  are  "structured"  (in  the  sense  that  the  states  are  of 
quadratic  form)  a  novel  approach  of  expressing  this  non-linear  dynamics  in  state  space 
form  is  given.  This  should  provide  a  useful  way  to  implement  the  equations  in  a 
computer  simulation  program  and  possibly  for  future  application  of  non-linear 
analysis  and  synthesis  techniques.  Detailed  consideration  of  the  quadratic/bilinear 
type  of  dynamic  systems  is  given  in  [4]. 

This  report  also  considers  a  locally  linearised  state  space  model  that  lends  itself  to 
better  known  linear  techniques  of  the  modem  control  theory.  A  coupled  multi-input 
multi-output  (MIMO)  model  is  derived  suitable  for  both  the  application  of  the  modem 
control  techniques  as  well  as  the  classical  time-domain  and  frequency  domain 
techniques.  This  is  validated  by  comparing  the  model  with  others  previously  published 
and  through  simulation  of  a  decoupled  single-input  single-output  (SISO)  system. 


2.  Rigid  Body  Dynamics 

2.1  Notation  and  Convention 

Conventions  and  notations  for  vehicle  body  axes  systems  as  well  as  the  forces, 
moments  and  other  quantities  are  shown  in  Figure  2.1  and  defined  in  Table  2.1. 


Figure  2.1  Motion  variable  notations 
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The  variables  shown  in  Figure  2.1  are  defined  as: 
m  -  mass  of  a  vehicle. 
a  -  incidence  in  the  pitch  plane. 

(]- incidence  in  the  yaw  plane. 

A  -  incidence  plane  angle. 

a-  total  incidence,  such  that:  tan  a  =  tan  a  cos  A,  and  tan  ft  =  tan  crsin  A. 
T  -  thrust. 


Table  2.1:  Motion  variables 


Vehicle  Body  Axes  System 

Roll 

axis 

X 

Pitch 

axis 

Y 

Yaw 

axis 

Z 

Angular  rates 

V 

R 

r 

Component  of  vehicle  velocity  along  each  axis 

u 

V 

w 

Component  of  aerodynamic  forces  acting  on  vehicle  along 
each  axis 

X 

Y 

Z 

Moments  acting  on  vehicle  about  each  axis 

L 

M 

N 

Moments  of  inertia  about  each  axis 

Ixx 

lyy 

In 

Products  of  each  inertia 

Iyz 

Lx 

I* 

Longitudinal  and  lateral  accelerations 

ax 

fly 

az 

Euler  angles 

9 

¥ 

Gravity  along  each  axis 

£*  . 

%y 

& 

Vehicle  thrust  along  the  body  axis 

T 

aileron  deflection. 
jj  -  elevator  deflection, 
g-  rudder  deflection. 


Figure  2.2  defines  the  control  surface  convention.  Here  the  control  surfaces  are 
numbered  as  shown  and  the  deflections  ( S ,,  52,  S3,  S4 )  are  taken  to  be  positive  if 
clockwise,  looking  outwards  along  the  individual  hinge  axis.  Thus: 

Aileron  deflection:  £  =  ~(81  +  5 2  +  S3  +  8 4  ) ,  if  all  four  control  surfaces  are  active;  or 

4 

%  =  —(5j  +53) ,  or  £  =  —  (82  +  54)  if  only  two  surfaces  are  active.  Positive  control 
2  2 
defection  (4)  causes  negative  roll. 
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Elevator  deflection:  r}  =  ~(8,-8t).  Positive  control  deflection  (r?)  causes  negative 
pitch. 

Rudder  deflection:  £  =  1-(82  -84).  Positive  control  deflection  (0  causes  negative  yaw. 


Figure  2.2  Control  surfaces  seen  from  the  rear  of  a  missile 


2.2  Euler's  Equations  of  Motion 

The  six  equations  of  motion  for  a  body  with  six  degrees  of  freedom  may  be  written  as 

[1-3]: 

m(u  +  wq-vr)  =  X  +  T  +  gxm  (2, ^ 

m(v  +  ur -wp)  =  Y  +  gym  (2-2) 

m(w—uq  +  vp)  —  Z  +  gxm  (2- 


I«p-(Iyy-D<ir  +  IJr2-q2)-IJpq  +  r)  +  Ixy(rp-q)  =  L 

Iyyq-(I«-Urp  +  IJP2-r2)-IJ<ir  +  P)  +  I>*(Pq-i')  =  M 
iJ-(i„-iyy)P‘i+IJ‘i2-P2)-IJrP+i)+I°‘(qr'p)  =  N- 


(2.4) 

(2.5) 

(2.6) 


Here 


-  is  the  derivative  operator. 
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Equations  (2.1  to  2.3)  represent  the  force  equations  of  a  generalised  rigid  body  and 
describe  the  translational  motion  of  its  centre  of  gravity  (c.g)  since  the  origin  of  the 
vehicle  body  axes  is  assumed  to  be  co-located  with  the  body  c.g. 

Equations  (2.3  to  2.6)  represent  the  moment  equations  of  a  generalised  rigid  body  and 
describe  the  rotational  motion  about  the  body  axes  through  its  c.g. 


Separating  the  derivative  terms  and  after  some  algebraic  manipulation.  Equations  (2.1 
to  2.3)  may  be  written  in  a  vector  form  as: 


d_ 

dt 


u 

V 

= 

w 

0 

-1 

0 


1 

o 

r-H 

O 

uq 

~X  +  T 

~8x 

0  0  10 

-10  0  0 

nr 

vp 

vr 

wp 

wq 

+ 

Y 

Z 

+ 

Sy 

_Sz_ 

(2.7) 


where:  X  - 


m 


Note  that  the  states  (u.V/W.p^r)  appear  as  "quadratic  terms  (form)". 


Equations  (2.4  to  2.6)  can  be  written  as: 


V 

/ 

~  L 

q 

-W 

pq 

+ 

M 

r 

pr 

N 

q2 

qr 

r2 

(2.8) 


where:  matrices  [a]  and  M  are  given  by: 
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I XX 

I  xy  I  zx 

M- 

-4 

I  yy  ~  4 

-4 

I yz  4 

0 

4 

-4 

4 

(4-4) 

-4 

M- 

-4 

1 

p 

1 

a 

0 

4 

4 

4 

V  xx  ^  yy  ) 

4 

-4 

-4 

0 

Here  again,  the  states  (p,q,r)  appear  in  "quadratic  form". 
Equation  (2.8)  may  also  be  written  as: 


p 

=k'Ml 

V 

+k'l 

'  4 

pq 

M 

r 

pr 

N_ 

q2 

qr 

r2 

where  the  inverse  W  is  given  by  (see  Appendix  A. 3): 


mM 


yz~  xy  I yyl zx  ) 


if  yy  I  zz  lyz  )  if zz^ xy  +  44f)  if yz^ xy  ~  yy zx  / 

if  zz  I  xy  +  44  )  if xx^  zz  ~  4  )  144  +  44-  ) 

)  [ijyz+lxylzx)  V  J yy  "V) 


if  yz^  xy  f  yy  ^  zx 


and  A 


-fc 


xx  ^ yyl ’  z: 


■I  I 

XX  yz 


■  T  T  2  -  T  I  2  -21  I  I 

1  yy1  zx  1  zz1  xy  yz1  zx1  xy 


)• 


(2.9) 


(2.9a) 


The  selection  of  the  particular  order  of  the  terms  in  the  "quadratic-state  vectors' 

\uq  ur  vp  vr  wp  wqY  ot  Equation  (2.7)  and  \p2 pq  pr  q2qr  r  °f  Equation  (2.8)  is 
discussed  in  Appendix  A.l. 


Combining  Equations  (2.7)  and  (2.9),  we  obtain  the  full  6th  order  rigid  body  dynamics 
state  equations  as: 
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d 

'-?r 

>1 

;  W 

_ ] _ 

V/r 

-j- 

W  j  M  ' 

«Y 

g 

dt 

£ 

[d 

IN"' Ml 

w  Ur] 

0 

(2.10) 


where  [c]  = 


0 

0 

1 


0  010-1 
-1001  0 
0-1000 


=[u  V  wf, 
x[?=\p  q  rf, 
jc[2]  =  [uq  ur  vp  vr  wp  wqf , 

hi  r  2  2  2 1^ 

X\  =  [p  pq  pr  q  qr  r  J  , 
u[^=[x  +  T  Y  zf, 
u[?=[L  M  Nj , 

g  =  [gx  Sy  gJ  • 

Equation  (2.10)  may  be  written  in  a  compact  form  as: 
=  [f]*^  +  [g]w^+  g^ 


(2.11) 


where 


[fi¬ 

fe]- 


‘[c] 

j  M  ' 
1 

Jo] 

i  [M14 

"M 

1 

1 

!  ur  l 

:  is  the  6x12  (quadratic)  state  coefficient  matrix. 


the  6x6  coefficient  matrix. 
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=[xr1'1  |  ;c[i'1]r  =[M  v  w  |  p  q  rj :  is  the  6x1  linear-state  vector. 

x1’1  =  [x1/1  U[;*]r  =  [uq  ur  vp  vr  wp  wq\  p'  pq  pr  q'  qr  r  ]  :  is 

the  12x1  quadratic-state  vector. 

Ml'l  =[x  +  f  Y  Z  \  L  M  tv]'  :  is  6x1  a  vector  function  of  control 

inputs,  forces  and  moments. 

and  g[,|=y  o\  =[gx  gy  g[  0  0  0]T:  is  the  6x1  gravity  (or  disturbance) 
vector. 

Note  that  for  a  two-axis  symmetrical  airframe,  ly2  =  Izx  =  Ixy  -  0.  Hence,  in  this  case,  the 
equation  (2.9)  can  be  reduced  to: 


p 

~0  0 

pq 

~L~ 

q 

= 

0  /„  0 

pr 

+ 

M 

r 

- 1 

<£> 

o 

1  a 
is 

_ _ I 

qr. 

N 

where  /„ 


L 


As  a  result,  the  state  equation  for  the  system  now  becomes: 


(2.12) 


if 

>] 

w 

"[/] 

[»] 

W 

g 

d 

= 

+ 

+ 

dt 

if. 

[0] 

W 

.if1. 

[0] 

[/] 

0 

(2.13) 


where:  \h]  - 


0 

0 


0  L 


yy 

0 


0 

0 


,  * =  \pq  Pr  qrT ' 
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u11]  =  [l  m  ivf. 

Remarks: 


Equations  (2.11)  and  (2.13)  are  complete  non-linear  description  of  the  full  6-DOF 
autopilot  model.  In  fact,  these  equations  contain  quadratic  terms  in  states  and  will  be 
classed  as  the  quadratic  dynamic  model.  This  type  of  model  is  required  when  autopilot 
design  is  undertaken  for  a  missile  executing  high  g-  or  high  angle  of  attack 
manoeuvres,  and  (u,  v,  w,  p,  q,  r)  are  not  small. 

A  more  detailed  consideration  of  the  algebraic  structure  of  this  type  of  dynamic 
systems  is  given  in  [4], 

2.3  Linearised  model  for  a  two-axis  symmetrical  airframe 


It  is  assumed  that  X ,  Y ,  Z ,  L  ,  M  and  N  are  functions  of  u,  v,  w,  p,  q,  r,  q  and  g. 
Using  first  order  linearisation  about  the  nominal  values  uo,  Vo,  wo,  po,  qo,  To,  £0,  qo  and  Co, 
and  defining  the  aerodynamic  derivatives  as: 


X.. 


X, 


_8X  ~  _  dX  ~  _  dX  9 

-  ~r~ >  A v  ~  77 '  Aw  -  77 '  A P 
du  8v  ow 

_8X  ~  _8X  ~  _8X 

“  /  ~  ~  /  2jL  r  “  / 

d£,  9  dq  ^  8 ^ 


dX  ~  _8X_  ~  _8X 

/  xx  _  /  XX  r  / 

8p  q  8q  r  dr 


y  _dY  ~  _d?  y  _d¥  y  y  =dl  ~  =8? 

dw  dv  ow  dp  8q  dr 

y  _dY_  ~  _df  ~  _af 

dt'Ir>~  dq'  d£' 


_8Z_  ~  _dZ  ~  _dZ 

"a.  '  A  "  a.  '  •>  ;  a,  ' 

dw  dv  dw 

_dZ  ~  _dZ  ~  _dZ 

_d£'  9~  dq '  C  ~8C' 


z. 


dZ  ~  _  dZ  ~  _  dZ 

_  /  Xj  / 

dp  9  dq  dr 


r  -®L  r  ~^L  7  ~®L  7  -®L  7  7 

U~du'  v  ~  dv  '  w~  dw'  P~dp'  9  ~  dq  '  r~  8r' 
~  _8L_  r  _  dT  ~  _8L_ 
dt’  '"dr,'  <-dt' 
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~  dM  ~  dM  n  dM  ~  dM  ~  dM  n  dM 

Mu  = - ,  Mv  =  ,  Mw  =— — ,  M  =— —  /  Mq  =  Mr  - 

du  dv  dw  dp  dq  dr 

~  5M  ~  oM  ~  5M 

~  SN  a  8N  ~  av  ~  aV  ~  dN  ~  _8N 

N" =  aT'  ^  =  av '  Nw=lw' Np  =!>p' N<  =  a? '  ^  "a7' 

~  a#  ~  a#  ~  aw 

<~°e' 

The  six  equations  of  motion  of  an  airframe  (using  equation  (2.13))  can  thus  be  written 


J  ti  =  r0dv  +  vQAr  -  q0Aw  -  w0Aq 

+  (xuAu  +  XvAv  +  XwAw  +  X  pAp  +  XgAq  +  XrAr  +  X.Alq  +  X^Ar\  +  XQAq) 
+  AT  +  Agx 

(2. 

Av  =  /J0dw  +  w0Ap-r0Au-u0Ar 

+  (YuAu  +  YvAv  +  YwAw  +  YpAp  +  YqAq  +  YrAr  +  Y%A^  +  Y^Av\  +  f?^)+ 


(2.14a) 


(2.14b) 


Aw  =  q0Au  +  u0Aq  -  p0Av  -  v^Ap 

+  (zuAu  +  ZvAv  +  ZwAw  +  ZpAp  +  ZqAq  +  ZrAr  +  Z \A^  +  Z^Ar]  +  ZqAq)+  Agz 


Ap  ^I^oAr  +  r0Aq) 

+  (ZuAu  +  LvAv  +  LwAw  +  LpAp  +  LqAq  +  LrAr  +  +  L^Ar\  +  LqAq) 


(2.14c) 


(2.14d) 


Aq  =  Iyy(r0Ap  +  p0Ar) 

+  {MuAu  +  MvAv  +  MwAw  +  MpAp  +  MqAq  +  MrAr  +  M^A%  +  M^Ar]  +  MqAq) 


(2.14e) 


Ar  =  Izz{p0Aq  +  q0Ap) 

+  (NuAu  +  NvAv  +  NwAw  +  NpAp  +  NqAq  +  NrAr  +  N^Aq  +  N^Ax\  +  NqAq) 


(2.14f) 
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Equations  (2.14a  to  2.14f)  may  be  written  in  a  matrix  notation  as: 


Au 

xu~ 

(r0+Xv) 

(~q0+XJ 

(-™0+xv) 

(v0+Xr) 

dv 

(~r0  +YU) 

K  ^ 

(Po+YJ 

(wo  +Yp) 

y,m 

(~u0+Yr) 

Aw 

(<lo+  4  ) 

(~Po  +  4J 

zw 

(~V0  +Zp) 

(uo  +Zq) 

4  ^ 

Ap 

4 

4 

4 

4 

(7.^+lj 

(T* xPo  +4 ) 

Aq 

Mu 

4 

(7y°+^P) 

(lyyPo  +Mr) 

Ar 

4 

4 

4 

(L  q0+xp) 

(IBPo+Nq) 

4  j. 

X,  Xn  Xr 


4 


L  4 


Mf  Mn  Mg 
Nc  N„  4 


~AT  +  Agx~ 

fAf 

Agy 

Ap 

+ 

ASz 

0 

Ac 

0 

0 

Au 

Av 

Aw 

Ap 

Aq 

Ar 


Note  that  Equation  (2.15)  is  represented  in  a  state-space  form  as: 

^-Ax  =  [F,]Ax  +  [Gi]Au1  +  Aw ’ 
at 


where  Ax  = 


Au 

Av 

Aw 

Ap 

Aq 


,  Au j  = 


At 

Ap 


Ar 


Aw, 


AT  +  Agx 

AS>- 

dg; 

0 

0 

0 


(2.16) 
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{r0+Xv) 

-q0+Xw 

-W0+X9 

v0  +Xr 

-r0  +yu 

4  _ 

Po  +  4 

w0  +  4 

-u0+Yr 

<7o  +  ^ u 

~Po+Zv 

4 

-v0  +  4 

Uq  +  4 

4  _ 

4 

4 

4 

4  ^ 

I XX  r0  4 

^  xxtfo  4 

Mu 

Mv 

Mw 

Vo  +&„ 

*4 

Vo  +Mr 

K 

K 

Nw 

Vo  +ftP  Vo +4 

4 

T  =f0+  AT  ,  gx  =  gl0  +  Agx,  gy  =  gy0  +  Agy,  gz  =  +  Agz 


2.4  Incorporation  of  accelerometer  and  gyro  measurement  model 

Generally,  not  all  state  variables  in  the  state  equation  are  accessible  or  measurable.  The 
common  measurement  variables,  in  most  missiles  or  airplanes,  are  the  angular  rate 
components  (roll  rate,  p,  pitch  rate,  q,  and  yaw  rate,  r)  and  the  acceleration  components 

(ttx,  tty,  ttz). 

Assuming  that  the  gyros  provide  ideal  readings  of  the  angular  rates,  we  get: 


Pm  =P 
qm=<i 


=  r 


(2.17a) 

(2.17b) 

.(2.17c) 


where  pm,  qm  and  rm  are  the  measured  body  rates.  Normally,  errors  due  to  drifts  and 
noise  are  included.  These  appear  as  additional  additive  terms  in  equations  (2.17a)  to 

(2.17c).  , 

In  contrast  to  the  readings  of  the  angular  rate  components,  the  readings  of  the 
acceleration  components  are  dependent  on  the  location  of  the  accelerometers,  w.r.t.  the 

c.g.  of  the  body. 
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The  acceleration  components  measured  at  point  O  (where  O  is  at  a  distance  of  dx,  dy 
and  dz  from  the  central  point  of  gravity,  c.g.,  along  x-,  y-  and  z-axis,  respectively),  may 
be  written  as: 

ax  -u  +  qw-rv  -  dx(q2  +  r2)  +  d  y(pq  ~  r)  +  d  z{pr  +  q)  (2.18a) 

av  =v  +  ru  —  pw  +  dx(pq  +  r)  —  dy(p~  +r2)  +  d2(qr  —  p)  (2.18b) 


az  =  w  +  pv-qu  +  dx(pr-q)  +  dy(qr  +  p)-dz(p 2  +q2) 


(2.18c) 


If  the  accelerometers  are  mounted  along  the  x-axis  (ie.  dy  -  dz  -  0)  which  is  usually  the 
case,  then  equations  (2.18a-c)  reduce  to: 


ax  =  u  + qw-rv-dx(q2  +  r2 )  =  X  +  T  +  gx  -dx(q~  +r2 ) 
av  =v  +  ru-pw  +  dx(  pq  +  r)  =  Y  +  gy  +  dx(  pq  +  r ) 
a2  =  w  +  pv-qu  +  dx(  pr-q)  =  Z  +g,  +dx(  pr-q) 


(2.19a) 

(2.19b) 

(2.19c) 


Note  that  the  right  hand  side  of  Equations  (2.19a)  to  (2.19c)  come  directly  from 
Equations  (2.1)  to  (2.3). 

Linearising  Equations  (2.19a)  to  (2.19c),  and  using  the  relationship  (2.15)  gives  us: 

Ay(t )  =  [ H ,  ]Ax(  t)  +  [j,  ]Au(  t)  +  Av(  t)  (2.20) 

where:  Ay(t)  =  [dp  A q  Ar  Aax  Aay  Aa  Jr :  is  the  output  vector, 

iT 


[H.h 


~0  0  0  X 

u 

Y  +N  d 
u  u  ' 

Z  -M  d 
u  u 

0  0  0  X 

V 

Y  +N  d 

V  V 

Z  -M  d 

V  V 

0  0  0  X 

w 

Y  +N  d 
w  w 

Z  -M  d 
w  w 

10  0  X 

p 

Y  +(q  +1  q  +N  )d, 
p  0  zz  0  p 

Z  +(r  -I  r  -M  )d 

P  0  yy  0  p 

0  1  0  X  -2q  d 
q  0 

Y  +(p  +1  p  +N  )dx 
q  0  zz  0  q 

Z  -M  d 
q  q  ' 

0  0  1  X  -2r  d 
r  0 

Y  +N  d 
r  r  ’ 

Z  +(p  -I  p  -M  )d 
r  0  yy  0  r 

is 


the  state  output  matrix. 
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0 
0 
0 
X 

n 

Y+Nd , 

Z  -Mdx 

rj  0  x 

measurement  matrix, 


0 

0 

0 

Y(+N(dx 
Z4  - Mtd , 


0 

0 

0 

X 

Y+Nd 

?  f  x 

Z  -M  d 

i  s  x 


:  is  the  matrix  related  to  inputs  in  the 


Aul(t)  =  [A%  At]  AgJ  :  is  the  input  vector,  and 


Av,(t) 


A  rT 

0  0  0  —  +  Agx  Agy  Agt 
m 


:  is  the  disturbance  vector. 


2.5  Linearised  model  of  the  airframe  including  fin  servos 

Assuming  that  the  servo  dynamics  for  the  aileron,  elevator  and  rudder  can  be 
described  adequately  by  a  second  order  lag  as: 


AS, 

"  52 

+ 

+  1 

Ar\ 

Md 

52 

+ 

+  1 

© 

iil 

“in 

*.c 

“  52 

2 

+ 

+  1 

(2.21a) 

(2.21b) 


(2.21c) 


where  A£d,  Arjd  and  A£d  are  the  demand  aileron,  elevator  and  rudder  deflection, 
respectively. 

ks4,  ksr,,  and  k5(  are  the  servo  gain  for  the  aileron,  elevator  and  rudder, 
respectively. 
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fisn'  aRd  are  the  damping  factor  for  the  aileron,  elevator  and  rudder, 
respectively. 

cosf,  cosn,  and  cosr  are  the  natural  frequency  for  the  aileron,  elevator  and  rudder, 
respectively. 

Equations  (2.21a)  to  (2.21c)  can  be  converted  into  differential  equations  as  follows: 


M  =  ' 

(2.22a) 

Afj  =  -at  /  Ar\  -  2  fi  nG)  nAf]  +  k ^cd sn'  Ar/ d ,  (2.22b) 

A£  =  -a /AC  -  2u<co,(AC  +  ks(coJ  AC d .  (2.22c) 

Hence,  the  state-space  model  for  the  autopilot  of  a  missile  including  the  servos  and 
airframe  is: 

Ax2(t)  =  \A2\Ax2(t)  +  \B2\Au2(t)  +  Aw2(t) ,  (2.23) 

where 


Ax 2  (t)  =  [Au  Av  Aw  Ap  Aq  Ar  At,  Ar\  At,  At,  Ar\  At^  , 

Note  that  the  aileron,  elevator  and  rudder  deflection  now  become  state  variables. 
Hence,  the  dimension  of  the  state  vector  is  increased  to  [12  x  1], 


Au2(0  =  [^d  Ant  zK ;J, 


The  inputs  are  now  the  demanded  aileron,  elevator  and  rudder  deflection. 
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- 

0 

0 

0 

0 

0 

0 

[F,] 

fo] 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-®sn2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(Here  matrices  Fi  and  Gi  are  the  same  as  in  Equation  (2.16)). 

0  0 

0  0 

0  0 

0  0 

0  0 

0  0 

0  0 

0  0 

0  0 

0  0 

0 

o  k(0)s(2 

The  output  vector  (or  the  measurement  equation)  is  given  by: 

4 y2(t)  =  H2Ax2(t)  +  Av2(t),  (2.24) 

where  Ay2(t)  -  \/ipm  Aqm  Arm  AaXm  Aaym  Aa2m\ , 


w= 


0 

0 

0 

0 

0 

0 

0 

0 

0 

k  ,o  / 

s 4 

0 

0 


16 


DSTO-TN-0449 


(Note  that  gyro  drift  and  noise  and  the  accelerometer  bias  may  be  added  to  the  right 
hand  side  of  Equation  (2.24)). 


0  0  0 
0  0  0 
0  0  0 
0  0  0  ' 

0  0  0 
0  0  0_ 

(Here  matrices  Hi  and  }i  are  the  same  as  those  in  Equation  (2.20)). 

I-  -!  T 

Ar r 

Av2  (t)=  000000000  —  +  Agx  Ag  Ag2  . 

L  m 

Figure  2.3  shows  the  block  diagram  of  an  open-loop  autopilot  which  contains  the  fin 
servos  and  airframe. 


r 

i 


Demand  T?tn 


Figure  2.3  A  block  diagram  of  an  open-loop  autopilot  system. 
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2.6  Lateral  Auto-pilot  design 

For  the  case  of  small  perturbation,  we  may  assume  that  (u0/  v0,  w0/  po,  qo,  r0)  are 
identically  zero.  In  this  case,  the  airframe  model  decouples  into  two  lateral  dynamics 
(pitch  and  yaw)  and  one  roll  dynamics.  We  will  consider  the  lateral  autopilot  dynamics 
to  validate  the  model  derived  in  this  report. 

Figure  2.4  shows  the  block  diagram  of  a  closed-loop  autopilot  system. 


Figure  2.4  A  block  diagram  of  a  closed-loop  auto-pilot  system. 


Ignoring  the  instrument  (gyro,  accelerometer)  dynamics,  the  measured  roll,  pitch  and 
yaw  angular  rates  (the  gyro  outputs)  can  be  expressed  as  inputs  to  the  gyros  multiplied 
the  gyro  gains,  Kgr,  K&  and  K&,  respectively.  Similarly,  the  measured  longitudinal 
acceleration,  aXl  and  lateral  accelerations,  ay  and  aZr  are  inputs  to  the  accelerometers 
multiplied  accelerometer  gains,  Kax,  Kiy  and  Ka2,  respectively.  The  accelerometer  gains 
affect  the  steady  state  response  and  may  be  set  to  1  for  transient  tests.  Rescaling 
accelerometer  gains,  after  selecting  gyro  gains,  allows  a  unity  gain  autopilot  to  be 
designed. 
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The  reference  signals,  generally  used  for  testing  the  transient  time  response  of  the 
autopilot,  are  the  desired  accelerations  in  yaw  direction,  %<*,  the  pitch  direction,  azn,  and 
roll  rate,  pd.  The  reference  roll  rate  is  kept  at  zero  to  assess  the  missile  dynamics  in  roll 
against  spurious  disturbances.  Hence,  the  reference  vector,  Ar,  can  be  written  as: 

Ar  =  [Apd  Aazd  Aayd]T , 

For  a  case  of  lateral  directional  control,  the  control  input  signal  for  the  fin  servos  can  be 
derived  as  follows: 


Arjd  =  Aazd  - KazAaz  - KgqAq  (2.25a) 

Agd  =  Aayd  -  KayAay  -  KgrAr  (2.25b) 

For  sake  of  simplicity,  pd  is  set  to  zero  since  this  case  only  considers  the  lateral 
directional  control.  As  a  result,  the  control  input  vector  can  be  written  as: 

Au2(t)  =  Ar(t)-KAy,(t),  (2.26) 

where  K  is  the  feedback  matrix  as  follows: 


~0 

0 

0 

0 

0 

0  " 

K  = 

0 

K 

SR 

0 

0 

0 

Ka 

0 

0 

K 

gr 

0 

K 

ay 

0 

3.  Verification  of  the  Developed  Model 


In  order  to  verify  the  developed  model,  the  state-space  model  (Equations  (2.23,  2.24 
and  2.26))  was  converted  into  transfer-function  form  using  Matlab  symbolic  toolbox 
(see  Appendix  B  and  Appendix  C.l)  for  comparison  with  the  results  already  published 

[4]- 


Consider  the  following  derivatives  and  variables  to  be  non-zero: 

Lu,  Lp,  Lr ,  M q ,  Mw ,  M n ,  N r,  Nv,  N ^ ,  Xu,  X p,  X ^ ,Yr,Yv,Y^,Zw,Zri,ZtI,u0,dx, 


we  obtain  the  transfer-function  between  the  roll  rate  and  the  aileron  deflection  as: 


p(s)  _  +  ( Lu  +  X ^  -  X UL^  ) 

J(7)~  S2-(LP+  Xu  )s  +  (XuLp  -  XpLu ) ' 


(2.27) 
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The  transfer-function  between  the  pitch  rate  and  the  elevator  deflection  as: 

q(s) _ M qs  —  MqZw  +  MwZn _  (2.28) 

^)~  s2  ~(Mq+ZJs  +  MqZw--MwZq-Mwu0  ' 

The  transfer-function  between  the  yaw  rate  and  the  rudder  deflection  as: 

r(s) _ +  NVY^  - N(YV~ _ ,  (2.29) 

~£<7)  ~  s2  -(Nr+YJs  +  Nvu0  +  NrYv  -  NJr  ' 

The  transfer-function  between  longitudinal  acceleration  and  the  aileron  deflection  as: 

ax(s)  ( XUX^  +  )s  -  XuX^Lp  +  X  pX^Lu  ^30) 

%(s)  s2  ~(XU+Lp)s  +  XuLp  -  XpLu 

The  transfer-function  between  lateral  acceleration  ay  and  the  rudder  deflection  as: 

ay(s)  (t  +N(dx)s2  +  (NvY(dx - N( Yvdx  +  NYr- NrYt l )s  +  NJcu0  - N( Yvu0 
£(s 7  ~  ~~  s2  - (Nr  +YJs  +  Nvu0  +NrYv~ 


And  the  transfer-function  between  lateral  acceleration  az  and  the  elevator  deflection  as: 

ajs)  (Zn  +MtJdx)s2  +(MwZ1ldx-MrlZwdx  +  MJJZq  -MqZr})s  +  Mr;Zwu0  -MwZvu0 
ij(s)  ~  s2  ~(Mq  +Zw)s-Mwu0  -MwZq  +MqZw 


Furthermore,  if  it  is  assumed  that  Xu  =  X(  =Yr  —  Lu  -  0 ,  the  transfer-function  between 
the  roll  rate  and  the  aileron  deflection  may  be  simplified  to: 

P( s)  _  h-  (2.33) 

Z(s)  s-Lp' 

The  transfer-function  between  the  yaw  rate  and  the  rudder  deflection  as: 
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r(s)  Nes  +  NJ(-N(Yv 

C(s)  s2  ~(Nr  +Yv)s  +  Nvu0  +  NrYv  ' 


(2.34) 


And  the  transfer-function  between  the  lateral  acceleration  in  yaw  axis,  ay,  and  the 
rudder  deflection  measured  at  the  c.g.  can  be  rewritten  as: 

a/s)  Y€s2-NrYcs  +  NvYcu0-NYvu0 
C(s)  s2  ~(Nr  +Yv)s  +  Nvu0  +  NrYv 

Equations  (2.33)  to  (2.35)  are  identical  to  those  of  the  transfer-functions  presented  in 
[P.Gamell  and  D.J.East  [Equations  (4.6-6),  (4.6-8)  and  (4.6-7)]. 

The  state-space  model  was  used  for  simulation  of  open-loop  and  closed-loop  responses 
for  a  typical  missile,  using  the  same  values  as  those  used  in  [4]  (see  Appendices  C  and 

D). 

Figures  2.5  and  2.6  show  the  lateral  accelerations  of  the  missile  due  to  a  step  input  to 
the  rudder  and  elevator,  respectively,  for  an  open  loop  simulation.  As  can  be  seen  from 
these  figures,  there  are  large  steady  state  errors.  However,  the  steady  state  errors  can 
be  reduced  with  a  feedback  loop  as  can  been  seen  in  Figures  2.7  and  2.8.  These 
simulation  results  are  similar  to  the  results  presented  in  [4]. 


Figure  2.5  Open  loop  simulation:  Lateral  autopilot  response  to  a  step  demand  acceleration 
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Time  [s] 

Figure  2.8  Closed  loop  simulation:  Lateral  autopilot  response  to  a  step  demand  acceleration. 

4.  Conclusions 

Both  the  non-linear  and  linearised  autopilot  models  have  been  derived  in  this  report. 
The  state-space  model  of  a  missile  autopilot  was  validated  by  comparing  the  model 
with  the  other  published  results,  and  through  both  open  and  closed-loop  systems 
simulation.  The  non-linear  dynamics  model  presented  as  structural  quadratic  algebraic 
system  is  novel  and  will  be  used  for  developed  non-linear  control  techniques  suitable 
for  missile  systems  high  g-  manoeuvres  and  operating  of  a  range  of  aerodynamics 
conditions.  The  models  developed  in  this  report  are  useful  for  further  research  on 
precision  optimum  guidance  and  control.  It  is  hoped  that  the  higher  order  model  with 
motion  and  inertial  coupling  will  provide  more  accurate  representation  of  missile 
autopilot  dynamics  and  should  be  used  for  adaptive  and  integrated  guidance  and 
control  of  agile  missiles. 
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Appendix  A:  State  Equation 

A.l.  Quadratic  State  Vector 

Since  the  state  equations  (Equations  (2.7)  and  (2.9))  include  the  quadratic  terms,  which 
are  comprised  of  all  the  combinations  of  state  variables,  we  define  two  separate  state 
vectors,  one  linear  state  vector,  and  one  quadratic  state  vector. 

Let  define  the  linear-state  vector  as: 


II 

rR" 

e 

-a 

(Al-1) 

x[‘]  =\p  q  rj , 

(A1.2) 

x^=[x^  \  x^J  -\u  v  wj  p  q  r\ , 

(Al-3) 

We  shall  consider  the  quadratic-state  vector  x^ corresponding  to  x^.  The  quadratic- 
state  vector  will  be  defined  as  a  vector  whose  elements  are  components  of  a 
homogeneous  quadratic  polynomial  of  these  states  taken  in  the  same  lexicographic 
order.  That  is,  the  quadratic-state  vector  may  be  written  as: 

x ^  =  \u2uv  uw  up  uq  ur\v2vw  vp  vq  vr\  w2wp  wq  wr\p2 pq  pr\q2 qr  \r  T  (A1.4) 

'  III 

There  are  21  terms  in  this  quadratic  state  vector.  Note  that  the  dimension  of  the 

n\n  + 1 ) 

quadratic  state  vector  is  v  — — -  when  the  dimension  of  the  linear-state  vector  is  n. 

This  type  of  representation  has  been  used  by  other  authors  [3,  4]  when  describing  high 
order  state  combinations  of  dynamical  systems.  In  the  rigid  body  dynamic  equations 
(2.7)  and  (2.9),  coefficients  of  a  number  of  these  terms  are  zero.  For  the  sake  of 
simplicity  (to  avoid  setting  large  number  elements  in  the  matrices  to  zero),  only  those 
quadratic  states  that  are  associated  with  non-zero  terms  are  retained.  That  is,  the 

quadratic-state  vector  x^and  its  partitioned  form  may  be  written  as: 
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\x[2]]  4  \ 


uq  ur 


vp  vr 


wp  wq  \p2 

i 


5) 

A.2.  Linearisation  of  the  Quadratic  V ector 


pq  pr  q2 


qr  r2  (Al. 


Given  the  quadratic-state  vector  (A1.4),  the  first-order  locally  linearised  vector  is  given 
by: 


~(u0Aq  +  q0Au ) 

~<lo 

0 

0 

0 

u0 

0  " 

( u0Ar  +  r0Au ) 

ro 

0 

0 

0 

0 

u0 

(v0Ap  +  p0Av) 

0 

Po 

0 

0 

0 

(v0Ar  +  r0Av) 

0 

ro 

0 

0 

0 

vo 

(w0Ap  +  p0Aw) 

0 

0 

Po 

w0 

0 

0 

(w0Aq  +  q0Aw) 

0 

0 

Po 

0 

Wq 

0 

2p0Ap 

0 

0 

0 

2p0 

0 

0 

{p0Aq  +  q0Ap) 

0 

0 

0 

Po 

Po 

0 

( Po Ar  +  ro^P) 

0 

0 

0 

ro 

0 

Po 

2q0Aq 

0 

0 

0 

0 

2p0 

0 

(q0Ar  +  r0Aq) 

0 

0 

0 

0 

r0 

Po 

2q0Aq 

_0 

0 

0 

0 

0 

2f0_ 

Au 

Av 

Aw 

Ap 

Aq 
Ar  _ 

(A2.1) 


Note  that  (uo,  vo,  wo,  po,  qo,  ro)  are  local  operating  states. 
We  shall  write  this  relationship  in  a  compact  form  as: 

*l?]=[XoW] 

The  matrix  [X0  ]  is  defined  via  the  equation  (A2.1). 


(A2.2) 


A.3.  Calculation  of  Inverse  Matrix 

Given  a  matrix  \g\  defined  as: 


A 

-F 

-E 


-F 

B 

-D 


-E 

-D 

C 
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Its  inverse  [c]  1  is  given  by: 

\{bC-D2)  (< CF  +  DE )  (. DF  +  BE ) 
[g]-;  =  L  (cf+de)  (ac-e2)  (AD  +  EF)  , 

A  (DF  +  BE)  (. AD  +  EF )  (AB-F2) 
where  A  =  (aBC  -  AD2  -BE2  -CF2  -  2DEF ). 
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Appendix  B:  From  State-Space  Form  into  Transfer- 


function 

Consider  a  system  given  as: 

x(t)  =  Ax(t)  +  Bu(t) ,  (B-1) 

and  y(t)  =  Cx(t)  +  Du(t) .  (B-2) 

Laplace  transformation  of  Equations  (B.l)  and  (B.2)  yields: 

sX(s)  =  AX(s)  +  BU(s) ,  (B-3) 

and  Y  (s)  =  CX(s)  +  DU(s) .  (B-4) 

Equation  (B.3)  can  be  rearranged  as: 

X(s)  =  (sI-A)-'BU(s).  (B-5) 

Substituting  (B.5)  into  (B.4),  we  obtain: 

Y(s )  =  [C(sl  -  A)'1  B  +  D  ]U  (s)  =  H(s)U(s).  (B.6) 

Hence,  the  transfer-function  of  the  system  is: 

H(s)  =  ^  =  C(sl  -  A)'1  B  +  D  .  (B-7) 
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Appendix  C:  Matlab  Codes 

C.l.  Converting  the  state-space  model  into  the  transfer-  functions 

The  following  m.file  converts  the  state-space  model  into  transfer-functions 


clear  all; 

%  Define  the  deravatives  and  parameters 

syms  Xu  Xp  Xxi  f; 

syms  Yv  Yr  Yxi  Yzeta  f; 

syms  Zw  Zq  Zeta  f; 

syms  Lu  Lp  Lxi  f; 

syms  Mw  Mq  Meta  f ; 

syms  Nv  Nr  Nxi  Neta  Nzeta  f; 

syms  u  v  w  p  q  r  f; 

syms  xi  eta  zeta  xi_dot  eta_dot  zeta_dot  f; 

syms  xi_d  eta_d  zeta_d  f; 

syms  ks  ws  mus  f; 

syms  ax  ay  az  f; 

syms  dx  uO  f; 

syms  s ; 

%  Define  the  elements  of  the  matrix  A 

all=Xu;  al2=0;  al3=0;  al4=Xp;  al5=0;  al6=0; 
a2 1=0 ;  a22=Yv;  a23=0;  a24=0;  a25=0;  a26=-uO+Yr; 
a31=0;  a32=0 ;  a33=Zw;  a34=0;  a35=uO+Zq;  a36=0; 
a41=Lu;  a42=0;  a43=0;  a44=Lp;  a45=0;  a46=0; 
a51=0;  a52=0;  a53=Mw;  a54=0;  a55=Mq;  a56=0; 
a61=0 ;  a62=Nv;  a63=0;  a64=0;  a65=0;  a66=Nr; 

%  Define  the  elements  of  the  matrix  B 

bll=Xxi ;  bl2=0;  bl3=0; 
b21=Yxi ;  b22=0;  b23=Yzeta; 
b31=0 ;  b32=Zeta;  b33=0; 
b41=Lxi;  b42=0;  b43=0; 
b51=0 ;  b52=Meta;  b53=0; 
b61=Nxi ;  b62=Neta;  b63=Nzeta; 

%  Define  the  elements  of  the  matrix  C 

c 1 1  =  0 ;  Cl2  =  0 ;  Cl3  =  0;  cl4  =  l;  Cl5  =  0;  cl6  =  0; 
c21=0 ;  C22=0 ;  c23=0;  c24=0;  c25=l;  c26=0; 
c31=0 ;  c32=0 ;  c33=0;  c34=0;  c35=0;  c36=l; 
c41=Xu;  c42=0 ;  C43=0;  c44=Xp;  c45=0;  c46=0; 

C  5 1  =  0 ;  c52=Yv+Nv*dx;  c53  =  0;  c54  =  0;  C55  =  0;  c56=Yr+Nr*dx; 
c61  =  0 ;  c  6  2  =  0 ;  c63  =  Zw+Mw*dx;  c64  =  0;  c65  =  Zq-Mq*dx;  c66=Yr; 
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%  Define  the  elements  of  the  matrix  D 

dll=0 ;  dl2=0 ;  dl3=0; 
d21=0 ;  d22=0 ;  d23=0; 
d3 1=0 ;  d32=0 ;  d33=0; 
d41=Xxi;  d42=0 ;  d43=0; 

d51=Yxi+Nxi*dx;  d52=0;  d53=Yzeta+Nzeta*dx; 
d61=0 ;  d62=Zeta-Meta*dx;  d63=0; 

A= [all , al2 , al3 , al4 , al5 , al6 ; 
a21 , a2  2 ,  a23 , a24 , a25 , a26 ; 
a31 , a32 , a33 , a34 , a35 , a36 ; 
a41 , a42 , a43 , a44 , a45 , a46; 
a51 ,  a52 , a53,a54,a55, a56 ; 
a61,a62,a63,a64,a65,a66]  ; 

B= [bll , bl2 , bl3 ; 
b21,b22,b23; 
b3 1 , b32 , b33 ; 
b41 , b42 , b43 ; 
b51 , b52 , b53 ; 
b61,b62,b63]  ; 


C= [ell , cl2 , cl3 , Cl4 , cl5 , Cl6 ; 
c21 , c22 , c23 , c24 , c25 , c26 ; 
C31,c32,c33,c34,c35,c36; 
C41,c42,c43,c44,c45,c46; 
051,052,053,054,055,056; 
C61,c62,c63,c64,c65,c66]  ; 

D= [dll , dl2 , dl3 ; 
d21,d22,d23; 
d3 1 ,  d3  2  ,  d3  3  ; 
d4 1 ,  d4  2  ,  d4  3  ; 
d51 , d52 , d53 ; 
d61, d62 ,  d63]  ; 

U= [xi ; eta ; zeta] ; 

S=  [s, 0,0, 0,0,0; 

0,s,0,0,0,0; 

0,0,3, 0,0,0; 

0,0, 0,3, 0,0; 

0,0, 0,0, s,0; 

0, 0, 0, 0, 0,s] ; 

%  The  transfer-function  of  the  system  is: 

%  H(s)  =  ^  =  C(SI-A)-1B  +  D 
u(s) 
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%  Calculate  (SI -A)  1 

SIAinv=  inv(S-A); 

%  Calculate  the  transfer- function 

H=C*SIAinv*B  +D; 

%  Display  the  transfer- function  between  the  lateral 
%  acceleration,  ay,  and  the  rudder  deflection, 

H  (5 , 3  ) 


C.l.  Open-loop  and  closed-loop  simulation 

%********************************************************* 

%This  main  program  simulates  the  response  of  the  missile 
%for  a  step  input  at  the  elevator  and  rudder. 

clear  all; 

%  Call  the  M-file  called  values .m 

values ; 

%  Call  the  M-file  called  init.m 

init  ; 

%  Call  the  M-file  called  ssmodel .m 

ssmodel ; 
for  t  =  0 : ts  :  5 

REF=  [theta_d;az_d,-ay_d]  ; 

U=REF  -  K*Y; 

X_dot_prev=X_dot ; 

X_dot=A*X+B*U; 

X=  X+ (X_dot_prev  +  X_dot)/2*ts; 

Y=C*X; 

accy= [accy, Y(5, 1)  ]  ; 
accz=  [accz, Y (6,1)]  ; 
time  =  [time,t]; 
ay_d=50 ; 
az_d=  50; 
theta_d  =  0 ; 
end; 

plot (time, accl, time, acc3 ) ; 
xlabel('Time  [s] ' ) ; 

ylabel (' Lateral  acceleration,  a_y  [m/s^2]'); 

title  (  1  Closed  loop  response  for  a  step  input  of  50  m/s'A'2  at  the 

rudder 1 ) ; 

figure; 

plot (time, acc2 , time , acc3 )  ; 
xlabel ( 1  Time  [s]  ' )  ; 
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ylabel ( 'Lateral  acceleration,  a_z  [m/sA2]'); 

title  (' Closed  loop  response  for  a  step  input  of  50  m/s  2  at  the 
elevator' ) ; 


%****************  Values. m  ******************** 
%  Set  values  of  the  parameters 

Xu=0 . 0 ;  Xp=0 . 0 ;  Xxi=0.0;  Xeta=0.0;  Xzeta=0.0; 


******* 


** 


*  * 


Yv= - 3 ;  Yp=0 ;  Yr=0;  Yxi=0.0;  Yeta=0;  Yzeta=180; 

Lu=  0.0;  Lp=0 . 0 ;  Lxi=0.0; 

Mu=0.0;  Mv=0.0;  Mw=-1.0;  Mq=-3.0;  Meta=-500.0;  Mzeta=0.0; 

Nv=  1;  Np=0 . 0 ;  Nw=0.0;  Nr=-3;  Nxi=0.0;  Neta=0.0; 
Nzeta=-500 . 0 ; 


Zw= - 3 ;  Zq=Yr ;  Zeta  =-Yzeta; 


ks=0 . 0068 ;  ksz=ks 
mus=0 . 7 ;  mus  z=mus ; 
ws=180 ; 
dx  =0.5; 
u0=500 ; 


%  For  the  case  of  open-loop  simulation ,  Kgr-Kay-Kgq-Kaz=0 

Kgr=  30.75;  Kay  =  0.825; 

Kgq=-Kgr;  Kaz  =  Kay; 


%***********************  init.m  ************************* 
%Initialise  the  parameters 

u=0 ;  v=0 ;  w=0;  p=0;  q=0;  r=0; 

xi=0;  eta=0 ;  zeta=0; 

xi_dot=0;  eta_dot=0;  zeta_dot=0; 

X=  [u; v; w; p ; q ; r ; xi ; eta ; zeta ; xi_dot ; eta_dot ; zeta_dot ]  ; 


X_dot=X; 

Y= [ 0 ; 0 ; 0 ; 0 ; 0  ;  0 ]  ; 


i  =  l; 
ay_d=0 ; 
az_d=0; 
theta_d=0 ; 
xi_d(i) =0; 
eta_d(i) =0; 
zeta  d (i) =0 ; 


ts=0 . 001; 
accy=0 ; 
accz=0 ; 
time=0 ; 
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%  *********************  ssmodel .m  *********************** 
%  Define  the  state-space  model 

al_l=Xu;  al_2=0  ;  al__3  =  0  ;  al_4=Xp ;  al_5=0  ;  al_6=0 ; 
al_7=Xxi ; al_8=0 ; al_9=0 ; al_10=0 ; al_ll=0 ;al_12=0; 

a2_l=0 ; a2_2=Yv; a2_3=0 ;a2_4=0 ; a2_5=0 ; a2_6=-uO+Yr ; 
a2_7=Yxi ; a2_8=0 ; a2_9=Yzeta; a2_10=0 ; a2_ll=0 ; a2_12=0 

a3__l  =  0  ;  a3_2=0 ;  a3_3 =Zw; a3_4=0  ;  a3_5=uO+Zq;  a3_6  =  0  ; 
a3_7=0 ; a3_8=Zeta; a3_9=0 ;a3_10=0 ; a3_ll=0 ; a3_12=0 ; 

a4_l=Lu; a4_2=0 ; a4_3=0 ;a4_4=Lp; a4_5=0 ; a4_6=0 ; 
a4_7=Lxi  ;  a4__8  =  0 ;  a4_9=0  ;a4_10  =  0  ;  a4_ll=0  ;  a4_12  =  0  ; 

a5_l=0 ; a5_2=0 ; a5_3=Mw; a5_4=0 ; a5_5=Mq; a5_6=0 ; 
a5_7=0 ; a5_8=Meta; a5_9=0 ; a5_10=0 ; a5_ll=0 ; a5_12=0 ; 

a6_l=0 ; a6_2=Nv;a6_3=0 ; a6_4=0 ; a6_5=0 ;a6_6=Nr ; 

a6_7=Nxi ; a6_8=Neta; a6_9=Nzeta; a6_10=0 ; a6_ll=0 ; a6_12=0 ; 

a 7 _ 1=0 ; a 7 _ 2  =  0 ; a 7 _ 3  =  0 ; a7_4  =  0 ; a 7 _ 5  =  0 ; a 7 _ 6  =  0 ; 

a7_7=0 ; a7_8=0 ; a7_9=0 ; a7_10=l ; a7_ll=0 ; a7_12=0 ; 

a8_l=0 ; a8  2=0 ;a8 _ 3=0 ;a8 _ 4=0 ;a8  5=0 ;a8 _ 6=0; 

a8_7=0 ; a8_8=0 ; a8_9=0 ; a8_10=0 ; a8_ll=l ; a8_12=0 ; 

a9_l=0 ; a9_2=0 ;a9_3=0 ; a9_4=0 ; a9_5=0 ; a9_6=0 ; 
a9_7=0 ; a9_8=0 ; a9_9=0 ; a9_10=0 ; a9_ll=0 ; a9_12=l ; 

al0_l=0;al0_2=0;al0_3=0;al0_4=0;al0_5=0;al0_6=0 ; 
al0_7  =  -ws^2 ; al0_8  =  0 ; al0_9=0 ; al0_10  =  -2  *mus*ws ; al0_ll=0 ; 
al0_12=0 ; 

all_l=0 ; all_2=0 ; all_3=0 ;all_4=0 ; all_5=0 ; all_6=0 ; 
all_7  =  0;all_8  =  -wsA2 ; all_9=0 ; all_10  =  0 ; all_ll  =  -2*musz*ws  ; 
all_12=0 ; 

al2_l=0 ; al2_2=0 ; al2_3=0 ;al2_4=0 ; al2_5=0 ; al2_6=0 ; 
al2_7=0 ; al2_8=0 ; al2_9= -ws^2 ; al2_10=0 ; al2_ll=0 ; 
al2_12=-2*mus*ws ; 

bl_l=0;bl_2=0;bl_3=0 ; 
b2_l  =  0 ;b2_2  =  0 ; b2_3  =  0 ; 
b3_l=0;b3_2=0;b3_3=0; 
b4_l=0;b4_2=0;b4_3=0; 
b5_l=0;b5_2=0;b5_3=0; 
b6_l=0 ;b6_2=0 ; b6_3=0 ; 
b7_l=0;b7_2=0;b7_3=0; 
b8_l=0 ;b8_2=0 ; b8_3=0 ; 
b9  1=0 ;b9  2=0;b9_3=0; 
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bl0_l=-ks*ws^2;bl0_2=0 ;bl0_3  =  0  ; 
bll_l=0;bll_2=-ksz*ws^2 ;bll_3=0 ; 
bl2_l=0;bl2_2=0;bl2_3=-ks*wsA2; 

cl_l=0 ; cl_2=0 ; cl_3=0 ; cl_4=l ; cl_5=0 ; cl_6=0 ; 
cl_7=0 ; cl_8=0 ; cl_9=0 ; cl_10=0 ; cl_ll=0 ; cl_12=0 ; 

c2_l=0 ; c2_2=0 ; c2_3=0 ; C2_4=0 ; c2_5=l ; c2_6=0 ; 

C2_7=0 ; c2_8=0 ; c2_9=0 ; c2_10=0 ; c2_ll=0 ; c2_12=0 ; 

c3_l=0;c3_2=0 ; c3_3=0 ; C3_4=0 ; c3_5=0 ; c3_6=l; 
c3_7=0 ; c3_8=0 ; c3_9=0 ; C3_10=0 ; c3_ll=0 ; c3_12=0 ; 

C4_l=xu; c4_2=0 ; c4_3=0 ; c4_4=Xp ; c4_5=0 ; C4_6=0 ; 

c4_7=Xxi;c4_8=Xeta;c4_9=Xzeta;c4_10=0;c4_ll=0;c4_12=0; 

c5  1=0; c5_2=Yv+Nv*dx;C5_3=0 ; c5_4=0 ;c5_5=0 ; c5_6=Yr+Nr*dx; 

C5~7=Yxi+Nxi*dx;c5_8=Yeta+Neta*dx;c5_9=Yzeta+Nzeta*dx; 

C5~10=0 ; C5_ll=0 ; C5_12=0 ; 

C6  1=0 ; c6_2=0 ; c6_3=Zw-Mw*dx; c6_4=0 ; c6_5=Zq-Mq*dx; c6_6=0 ; 
c6_ 7=0 ; c6_8  =  Zeta-Meta*dx; c6_9  =  0 ; c6_10  =  0 ; c6_ll  =  0 ; c6_12=0 ; 

A= [al  1 , al_2 , al_3 , al_4 , al_5 , al_6 , al_7 , al_8 , al_9 , al_10 , al_ll , al_12 
a2  1 , a2_2 , a2_3 , a2_4 , a2_5 , a2_6 , a2_7 , a2_8 , a2_9 , a2_10 , a2_ll , a2_12 
a3  1 , a3_2 , a3_3 , a3_4 , a3_5 , a3_6 , a3_7 , a3_8 , a3_9 , a3_10 , a3_ll , a3_12 

a4  l,a4  2 , a4_3  , a4_4 , a4_5 , a4_6 , a4_7 , a4_8 , a4_9 , a4_10 , a4_ll, a4_12 

a5  1, a5_2 , a5_3 , a5_4 , a5_5 , a5_6 , a5_7 , a5_8 , a5_9 , a5_10 , a5_ll, a5_12 
a6  1 , a6_2 , a6_3 , a6_4 , a6_5 , a6_6 , a6_7 , a6_8 , a6_9 , a6_10 , a6_ll , a6_12 
a7  1, a7_2 , a7_3 , a7_4 , a7_5 , a7_6 , a7_7 , a7_8 , a7_9 , a7_10 , a7_ll , a7_12 
a8  1, a8_2 , a8_3 , a8_4 , a8_5 , a8_6 , a8_7 , a8_8 , a8_9 , a8_10 , a8_ll , a8_12 
a9_l^a9_2,a9_3,a9_4,a9_5,a9_6/a9_7/a9_8/a9_9,a9_10,a9_ll,a9_12; 

al0_l/al0_2,al0_3,al0_4^10_5/al0_6/al0_7,al0_8,al0_9/al0_10,al0_ll/ 
al0_12;  all_l/all_2/all_3/all_4,all_5/all_6/all_7,all_8,all_9,all_10,all_ll, 
all_12;  al2_l,al2_2/al2_3,al2_4/al2_5,al2_6,al2_7/al2_8/al2_9/al2_10/al2_ll/ 

al2_12]; 

B=  [bl_l,  bl_2  , bl_3  ; 
b2_l , b2_2 , b2_3 ; 
b3_l , b3_2 , b3_3 ; 
b4__l ,  b4_2  ,  b4_3  ; 
b5_l , b5_2 , b5_3 ; 
b6_l , b6_2 , b6_3 ; 
b7_l , b7_2 , b7_3 ; 
b8_l , b8_2 , b8_3 ; 
b9_l , b9_2 , b9_3 ; 
bl0_l , bl0_2 , bl0_3 ; 
bll_l , bll_2 , bll_3 ; 
bl2  1 , bl2_2 , bl2_3 ] ; 
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C= [ c 1 _ 1 , cl_2 , cl_3 , cl_4 , cl_5 , cl_6 

c2_l , c2_2 , c2_3 , c2_4 , c2_5 , c2_6 
C3_l , c3_2 , c3_3 , c3_4 , c3_5 , c3_6 
c4_l , c4_2 , C4_3 , c4_4 , c4_5 , c4_6 
c5_l , c5_2 , c5_3 , c5_4 , c5_5 , c5_6 
c6_l , c6_2 , c6_3 , c6_4 , c6_5 , c6_6 


D=  [0,0,0; 

0,0,0; 

0,0,0; 

0,0,0; 

0,0,0; 

0,0,0]  ; 

K  =  [0, 0, 0,  0, 0, 0; 

0 , Kgq ,0,0,0, Kaz ; 
0, 0, Kgr, 0, Kay, 0] ; 


cl_7 , cl_8 , cl_9 , cl_10  ,  cl_ll , cl_12 ; 
c2_7 , c2_8 , c2_9 , c2_10 , c2_ll , c2_12 ; 
C3_7 , c3_8 , c3_9 , C3_10 , c3_ll , c3_12 ; 
c4_7  ,  c4__8  ,  c4__9  ,  c4_10  ,  c4_ll ,  c4_12  ; 
c5_7 , c5_8 , c5_9 , c5_10 , c5_ll , c5_12 ; 
c6_7 , c6_8 , c6_9 , c6_10  ,  c6_ll , c6_12 ]  ; 
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Appendix  D:  Values  of  the  Non-Zero  Derivatives  and 
Parameters  used  in  the  simulation 


Yv  =  -3,  7,  =  180,  Nv  =  1,  Nr  =  -3,  Nc  =  -500, 

Zw  =  -3,  Z n  =  -180,  Mw  =  -l,  Mq  =  -3„  Mn  =  -500, 
ks4  =  0.0068,  kSJJ  =  0.0068,  ks(  =  0.0068, 

=0.7,  Ms,  =0.7,  =0.7, 

o)s4  =  180,  o)srj  =  180,  g)s4  =180, 
u0  =500,  dx  =  0.5 

Feedback  gains: 

Kgr=  30.75;  Kay  =  0.825; 

Kgq=-Kgr;  Kaz  =  Kay; 
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